################################################################################
#
#  REPLICATION SCRIPT
#
#  Title:   Disentangling the Three Facets of Mass Ideological Polarization:
#           A Network Approach across 78 Societies
#  Authors: Yufan Guo, Yilang Peng, Tian Yang
#  Journal: Public Opinion Quarterly
#
#  This script replicates ALL figures and tables in the main text and the
#  Supplementary Information (SI), in document order.
#
#  ─────────────────────────────────────────────────────────────────────────────
#  DOCUMENT MAP
#
#  MAIN TEXT
#    Figure 1  – Simulation: Correlation vs. Partial-Correlation Network
#    Figure 3  – World maps: global distribution of the three manifestations
#    Figure 4  – Scatter plots between the three manifestations
#    Table 2   – Top 10 societies in each manifestation
#    Figure 5  – OLS regression coefficient plot
#    Figure 6  – Scatter plots: HDI vs. each manifestation
#    Table 3   – OLS models predicting Disagreement
#
#  SUPPLEMENTARY INFORMATION
#    Section A – Full item list (Table A1) and country list (Table A2)
#    Section B – Robustness: Excluding governmental items (Figure B1)
#    Section C – Country-level variable descriptions (Tables C1–C2)
#    Section D – Robustness: Country-year as unit of analysis (Figure D1)
#    Section E – Measurement validation (Figures E1–E2)
#    Section F – Robustness: Within-operational alignment by density (F1–F3)
#    Section G – Robustness: Modernization by GDP per capita (Figures G1–G2)
#    Section H – Simulation detail (Figure H1)
#    Section I – Robustness: Weighted results (Figures I1–I2, Table I1)
#
#  ─────────────────────────────────────────────────────────────────────────────
#  REQUIRED DATA FILES  (place in the same folder as this script)
#
#  1. EVS_WVS_Joint_Rdata_v5_0.rdata
#       Joint EVS/WVS Wave 7 (2017–2022), version 5.0.0
#       Download: https://www.worldvaluessurvey.org/WVSEVSjoint2017.jsp
#
#  2. countrycode.xlsx  (provided in the replication archive)
#       Sheet 1 – maps WVS numeric country code (cntry_AN) to country names
#       Sheet 2 – country-level covariates: HDI, GDP (log), political stability,
#                 political parallelism, total population (log), press freedom,
#                 ethnic fractionalization, mean political interest (minterest),
#                 and V-Dem merge keys (country_text_id, year)
#
#  3. V-Dem-CY-Full+Others-v14.rds
#       V-Dem Country-Year full dataset, version 14
#       Download: https://v-dem.net/data/the-v-dem-dataset/
#
#  4. dimension_summary.csv  (provided in the replication archive)
#       Country-level scores for all three manifestations and their alternative
#       measures, pre-computed from the network estimation below.
#       Load this file to skip the computationally intensive loops (Sections
#       4–5) and proceed directly to figure replication (Section 6 onward).
#
#  ─────────────────────────────────────────────────────────────────────────────
#  REQUIRED R PACKAGES
#
#  install.packages(c(
#    "tidyverse", "psych", "bootnet", "qgraph", "igraph", "readxl",
#    "labelled", "agrmt", "questionr", "moments", "corrplot", "factoextra",
#    "cowplot", "ggrepel", "ggpubr", "sf", "cshapes", "viridis",
#    "fuzzyjoin", "texreg", "MASS", "scales"
#  ))
#
#  ─────────────────────────────────────────────────────────────────────────────
#  NOTE ON COMPUTATION TIME
#
#  The network estimation loops (Section 4) call estimateNetwork() for each
#  of the 78 countries.  Estimated run time: 10–30 minutes, depending on
#  hardware.  To reproduce figures without re-running, load the provided
#  dimension_summary.csv (Section 5) and skip Section 4.
#
################################################################################


################################################################################
# SECTION 1 — SETUP
################################################################################

# Set working directory to the folder that contains all data files.
# Uncomment and edit the line below for your local path:
# setwd("path/to/your/replication/folder")

library(tidyverse)    # Data wrangling and ggplot2
library(psych)        # Descriptive statistics (describe())
library(bootnet)      # Network estimation: estimateNetwork()
library(qgraph)       # Network visualization; centrality()
library(igraph)       # Graph metrics: edge_density()
library(readxl)       # Read Excel (.xlsx) files
library(labelled)     # Remove value labels from haven-imported data
library(agrmt)        # Tastle–Wierman consensus/disagreement: consensus()
library(questionr)    # Weighted frequency tables: wtd.table()
library(moments)      # Skewness and kurtosis for bimodality coefficient
library(corrplot)     # Correlation matrix visualization
library(factoextra)   # PCA visualization (SI Section E)
library(cowplot)      # Composite figures: ggdraw(), draw_plot_label()
library(ggrepel)      # Non-overlapping text labels in scatter plots
library(ggpubr)       # Multi-panel arrangement: ggarrange()
library(sf)           # Simple features (spatial data)
library(cshapes)      # World map shapefiles: cshp()
library(viridis)      # Perceptually uniform colour palettes: magma()
library(fuzzyjoin)    # Fuzzy/exact string joins for map data
library(texreg)       # Regression tables: screenreg()
library(MASS)         # Multivariate normal simulation: mvrnorm()
library(scales)       # Axis scale helpers: squish()


################################################################################
# SECTION 2 — DATA PREPARATION
#
# Load the joint EVS/WVS Wave 7 dataset, select the 34 operational ideology
# items plus the IDEOLOGY item (E033), recode response categories, replace
# invalid codes with NA, and define the analytic sample.
################################################################################

# 2a. Load the EVS/WVS joint dataset ──────────────────────────────────────────
load("EVS_WVS_Joint_Rdata_v5_0.rdata")
wvs <- `EVS_WVS_Joint_R_v5_0`
wvs <- labelled::remove_labels(wvs)   # strip haven value labels

# 2b. Load country code crosswalk (Sheet 1 of countrycode.xlsx) ───────────────
# Maps the numeric WVS country identifier (cntry_AN) to readable country names.
code <- read_excel("countrycode.xlsx", sheet = "Sheet1")

# 2c. Select and rename the 35 items used in the analysis ─────────────────────
#
# IDEOLOGY item (added to the 35-node network for symbolic-operational alignment)
#   E033  – Self-positioning on a 10-point left–right scale (1=left, 10=right)
#
# 34 OPERATIONAL IDEOLOGY ITEMS (used in both network types; see SI Table A1)
#
#  #  WVS code   Variable name    Item description
#  1  A124_02    diffrace         Neighbors: People of a Different Race (0/1)
#  2  A124_06    immigr           Neighbors: Immigrants/Foreign Workers (0/1)
#  3  B008       envir_econ       Protecting environment vs. Economic growth (1–3)
#  4  C001       jobscarce1       Jobs scarce: Men should have more right to a job (1–3, recoded)
#  5  C002       jobscarce2       Jobs scarce: Priority to nationals over immigrants (1–3, recoded)
#  6  D059       womanleader      Men make better political leaders than women (1–4)
#  7  D060       womenuniv        University more important for a boy than a girl (1–4)
#  8  D061       childsuffer      Children suffer when mother works full-time (1–4)
#  9  D078       womenbusin       Men make better business executives than women (1–4)
# 10  E035       incequal         Income equality (1=more equal, 10=less equal)
# 11  E036       ownership        Private vs state ownership of business (1–10)
# 12  E037       govresp          Government responsibility for citizens (1–10)
# 13  E039       competition      Competition: good or harmful (1–10)
# 14  E114       strongleader     Political system: Having a strong leader (1–4)
# 15  E115       experts          Political system: Having experts decide (1–4)
# 16  E117       dompol           Political system: Having democratic politics (1–4)
# 17  E224       taxrich          Democracy: Governments tax the rich (1–10)
# 18  E225       relig            Democracy: Religious authorities interpret laws (1–10)
# 19  E226       election         Democracy: People choose leaders in free elections (1–10)
# 20  E227       aidunempl        Democracy: People receive state aid for unemployment (1–10)
# 21  E228       army             Democracy: Army takes over when govt incompetent (1–10)
# 22  E229       civilright       Democracy: Civil rights protect liberty (1–10)
# 23  E233       womenright       Democracy: Women have same rights as men (1–10)
# 24  E233A      stateequal       The State makes people's incomes equal (1–10)
# 25  E233B      rule             People obey their rulers (1–10)
# 26  E235       impdem           Importance of democracy (1–10)
# 27  F118       homosex          Justifiable: Homosexuality (1–10)
# 28  F120       abortion         Justifiable: Abortion (1–10)
# 29  F121       divorce          Justifiable: Divorce (1–10)
# 30  F123       suicide          Justifiable: Suicide (1–10)
# 31  F144_02    deathpen         Justifiable: Death penalty (1–10)
# 32  G007_35_B  trustt1          Trust: People of another religion (1–4)
# 33  G007_36_B  trust2           Trust: People of another nationality (1–4)
# 34  G052       impact           Evaluate the impact of immigrants on the country (1–5)

wvs1 <- wvs %>%
  dplyr::select(
    wave, cntry_AN, year, gwght,    # survey metadata + weight variable
    E033,                           # IDEOLOGY item
    B008, C001, C002, D059, D060, D078, E035,
    E036, E037, E224, E225, E226, E227, E228, E229,
    E233, F118, F120, F121, F123,
    A124_02, A124_06, D061, E039, E114, E115, E117,
    E233A, E233B, E235, F144_02, G052, G007_35_B, G007_36_B
  ) %>%
  left_join(code, by = c("cntry_AN" = "countrycode")) %>%
  dplyr::select(-cntry_AN) %>%
  relocate(country, .after = wave) %>%
  dplyr::rename(
    ideology    = E033,      # left–right self-placement
    envir_econ  = B008,      # environment vs. economy
    jobscarce1  = C001,      # jobs scarce: men vs. women
    jobscarce2  = C002,      # jobs scarce: nationals vs. immigrants
    womanleader = D059,
    womenuniv   = D060,
    womenbusin  = D078,
    incequal    = E035,
    ownership   = E036,
    govresp     = E037,
    taxrich     = E224,
    relig       = E225,
    election    = E226,
    aidunempl   = E227,
    army        = E228,
    civilright  = E229,
    womenright  = E233,
    homosex     = F118,
    abortion    = F120,
    divorce     = F121,
    suicide     = F123,
    diffrace    = A124_02,
    immigr      = A124_06,
    childsuffer = D061,
    competition = E039,
    strongleader= E114,
    experts     = E115,
    dompol      = E117,
    stateequal  = E233A,
    rule        = E233B,
    impdem      = E235,
    deathpen    = F144_02,
    impact      = G052,
    trustt1     = G007_35_B,
    trust2      = G007_36_B
  ) %>%
  mutate(
    # C001/C002: original 3-category code has 1=men/nationals, 2=neither, 3=women/immigrants.
    # Swap 2 and 3 so that higher values consistently reflect more egalitarian attitudes.
    jobscarce1 = case_when(
      jobscarce1 == 1 ~ 1, jobscarce1 == 2 ~ 3,
      jobscarce1 == 3 ~ 2, TRUE ~ NA_real_),
    jobscarce2 = case_when(
      jobscarce2 == 1 ~ 1, jobscarce2 == 2 ~ 3,
      jobscarce2 == 3 ~ 2, TRUE ~ NA_real_),
    # B008: recode from 1=environment / 2=economy to binary 1/0
    envir_econ = case_when(
      envir_econ == 1 ~ 1, envir_econ == 2 ~ 0, TRUE ~ NA_real_),
    # E224–E229, E233, E233A–B, E235: "essential for democracy" items.
    # Code 0 is an invalid skip code in WVS; replace with NA.
    taxrich    = ifelse(taxrich    == 0, NA, taxrich),
    relig      = ifelse(relig      == 0, NA, relig),
    election   = ifelse(election   == 0, NA, election),
    aidunempl  = ifelse(aidunempl  == 0, NA, aidunempl),
    army       = ifelse(army       == 0, NA, army),
    civilright = ifelse(civilright == 0, NA, civilright),
    womenright = ifelse(womenright == 0, NA, womenright),
    stateequal = ifelse(stateequal == 0, NA, stateequal),
    rule       = ifelse(rule       == 0, NA, rule),
    impdem     = ifelse(impdem     == 0, NA, impdem)
  )

# Replace all remaining negative values (WVS "don't know", "not applicable",
# "no answer" codes: −1, −2, −3, −4, −5) with NA across all columns.
wvs1[wvs1 < 0] <- NA

# 2d. Define the analytic sample ──────────────────────────────────────────────
# Keep only respondents with a valid ideology self-placement (E033 ≠ NA).
# As noted in Appendix A, 12 societies did not measure self-placed ideology
wvs7 <- wvs1 %>% filter(!is.na(ideology))

# Reference vectors used throughout the analysis
countrylist <- as.character(unique(wvs7$country))  # 80 societies
# two (Tajikistan, Turkey) lacked one issue item — leaving 78 societies total.
nodelist    <- names(wvs7)[!names(wvs7) %in%
                 c("wave","country","year","gwght")]       # ideology + 34 items (35 nodes)
issuelist   <- nodelist[nodelist != "ideology"]            # 34 operational ideology items only


################################################################################
# SECTION 3 — HELPER FUNCTIONS
#
# One function per measure.  Each function takes a country name (n) and
# returns the computed measure for that country, wrapped in tryCatch() so
# that failures (e.g., too few observations) return NA without halting the loop.
################################################################################

# ── 3A. Disagreement (ordinal dispersion) ─────────────────────

disagreement_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)

    # Inner: item-level disagreement via agrmt::consensus()
    item_disagree <- function(item_name) {
      x <- dat[[item_name]]
      x <- x[!is.na(x)]
      if (length(x) < 10) return(NA_real_)
      freq_table <- as.vector(round(
        questionr::wtd.table(x, weights = rep(1, length(x))), 0))
      1 - agrmt::consensus(freq_table)
    }

    item_scores <- sapply(issuelist, item_disagree)
    mean(item_scores, na.rm = TRUE)
  }, error = function(e) NA_real_)
}

# Item-level disagreement (for SI Section A Table A1 correlations and
# Section E validation and issue-level comparisons)
disagreement_itemlevel_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)

    sapply(issuelist, function(item_name) {
      x <- dat[[item_name]]
      x <- x[!is.na(x)]
      if (length(x) < 10) return(NA_real_)
      freq_table <- as.vector(round(
        questionr::wtd.table(x, weights = rep(1, length(x))), 0))
      1 - agrmt::consensus(freq_table)
    })
  }, error = function(e) rep(NA_real_, length(issuelist)))
}

# ── 3B. Within-operational ideological alignment (sum of edge weights) ────────

within_op_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]

    g <- estimateNetwork(dat, default = "mgm")
    adj <- abs(g$graph)
    sum(adj[upper.tri(adj)])   # sum of edge weights
  }, error = function(e) NA_real_)
}

# ── 3C. Symbolic-operational ideological alignment (ideology InDegree) ────────

symbolic_op_fun <- function(n) {
  tryCatch({
    # Include both ideology and the 34 issue items (35 nodes total)
    dat <- subset(wvs7, country == n)[, c("ideology", issuelist), drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]

    g  <- estimateNetwork(dat, default = "mgm")
    ct <- centrality(g)
    # InDegree centrality vector: first element corresponds to the IDEOLOGY node
    ct[["InDegree"]][1]
  }, error = function(e) NA_real_)
}

# ── 3D. Alternative measure: Network density (SI Section F) ──────────────────

density_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]

    g     <- estimateNetwork(dat, default = "mgm")
    adj   <- abs(g$graph)
    graph <- igraph::graph_from_adjacency_matrix(adj, mode = "undirected",
                                                  weighted = TRUE)
    igraph::edge_density(graph)
  }, error = function(e) NA_real_)
}

# ── 3E. Alternative measure: Average pairwise Spearman correlation (SI Sec E) ─

avg_pairwise_cor_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    cor_mat <- abs(cor(dat, method = "spearman", use = "pairwise.complete.obs"))
    mean(cor_mat[upper.tri(cor_mat)])
  }, error = function(e) NA_real_)
}

# ── 3F. Alternative measure: Average |Spearman ρ| with ideology (SI Sec E) ───

avg_cor_ideology_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, c("ideology", issuelist), drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    ref  <- dat[["ideology"]]
    cors <- sapply(issuelist, function(v)
      abs(cor(ref, dat[[v]], method = "spearman", use = "pairwise.complete.obs")))
    mean(cors, na.rm = TRUE)
  }, error = function(e) NA_real_)
}

# ── 3G. Alternative measure: Bimodality coefficient (SI Section E) ───────────

bc_item_fun <- function(x) {
  x <- x[!is.na(x)]
  n <- length(x)
  if (n < 4) return(NA_real_)
  m3 <- moments::skewness(x)
  m4 <- moments::kurtosis(x) - 3          # excess kurtosis
  (m3^2 + 1) / (m4 + 3 * (n - 1)^2 / ((n - 2) * (n - 3)))
}

bimodality_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    mean(sapply(dat, bc_item_fun), na.rm = TRUE)
  }, error = function(e) NA_real_)
}

# ── 3H. Alternative measure: Mean variance across items (SI Section E) ────────

avg_var_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    mean(sapply(dat, var, na.rm = TRUE))
  }, error = function(e) NA_real_)
}


################################################################################
# SECTION 4 — RUN NETWORK ESTIMATION  [COMPUTATIONALLY INTENSIVE, ~10–30 min]
#
# This section computes all three manifestation scores and their alternative
# measures for all 78 societies.
################################################################################

cat("Computing Disagreement (ordinal dispersion)...\n")
disagreement_scores <- sapply(countrylist, disagreement_fun)

cat("Computing Within-operational ideological alignment (sum of edge weights)...\n")
within_op_scores <- sapply(countrylist, within_op_fun)

cat("Computing Symbolic-operational ideological alignment (ideology InDegree)...\n")
symbolic_op_scores <- sapply(countrylist, symbolic_op_fun)

cat("Computing alternative measures for validation...\n")
density_scores    <- sapply(countrylist, density_fun)
avg_pcor_scores   <- sapply(countrylist, avg_pairwise_cor_fun)
avg_ideo_scores   <- sapply(countrylist, avg_cor_ideology_fun)
bimodality_scores <- sapply(countrylist, bimodality_fun)
avg_var_scores    <- sapply(countrylist, avg_var_fun)

cat("Computing item-level disagreement for all countries...\n")
item_disagree_mat <- sapply(countrylist, disagreement_itemlevel_fun)
# Result: 34 × 78 matrix; rows = items (issuelist), columns = countries

# Combine into a single data frame
dimension_summary <- data.frame(
  country           = countrylist,
  disagreement      = disagreement_scores,
  symbolic_op_align = symbolic_op_scores,
  within_op_align   = within_op_scores,
  density           = density_scores,
  avg_pairwise_cor  = avg_pcor_scores,
  avg_cor_ideology  = avg_ideo_scores,
  bimodality_coef   = bimodality_scores,
  avg_variance      = avg_var_scores,
  stringsAsFactors  = FALSE
) %>%
  filter(!is.na(disagreement) & !is.na(symbolic_op_align) & !is.na(within_op_align))

# Save to CSV so downstream analyses can be reproduced without re-running
write.csv(dimension_summary, "dimension_summary.csv", row.names = FALSE)
cat("Saved dimension_summary.csv  (", nrow(dimension_summary), "societies)\n")

# Save item-level disagreement matrix
item_disagree_df <- as.data.frame(item_disagree_mat)
item_disagree_df$item <- issuelist
write.csv(item_disagree_df, "item_disagreement.csv", row.names = FALSE)
cat("Saved item_disagreement.csv\n")


################################################################################
# SECTION 5 — LOAD PRE-COMPUTED RESULTS
#
# Start here if you want to reproduce all figures without re-running Section 4.
# The file dimension_summary.csv is provided in the replication archive.
################################################################################

polar <- read.csv("dimension_summary.csv", stringsAsFactors = FALSE)
# Required columns:
#   country           – country name
#   disagreement      – Disagreement (ordinal dispersion)
#   symbolic_op_align – Symbolic-operational ideological alignment (ideology InDegree)
#   within_op_align   – Within-operational ideological alignment (sum of edge weights)
#   density           – Network density (alternative within-op measure; SI Sec F)
#   avg_pairwise_cor  – Avg pairwise Spearman ρ (alternative within-op; SI Sec E)
#   avg_cor_ideology  – Avg |Spearman ρ| with ideology (alternative symbolic-op; SI Sec E)
#   bimodality_coef   – Bimodality coefficient (alternative disagreement; SI Sec E)
#   avg_variance      – Mean variance across items (alternative disagreement; SI Sec E)

# Spearman correlations among the three main manifestations (reported in text)
cat("\n=== Inter-manifestation Spearman correlations (main text, Figure 4) ===\n")
cor.test(polar$disagreement, polar$symbolic_op_align, method = "spearman")
cor.test(polar$disagreement, polar$within_op_align,   method = "spearman")
cor.test(polar$symbolic_op_align, polar$within_op_align, method = "spearman")


################################################################################
# SECTION 6 — MAIN TEXT TABLE 2: Top 10 Societies in Each Manifestation
################################################################################

cat("\n=== Table 2: Top 10 societies – Disagreement ===\n")
top10_disagree <- polar %>%
  arrange(desc(disagreement)) %>%
  slice(1:10) %>%
  dplyr::select(country, disagreement)
print(top10_disagree)

cat("\n=== Table 2: Top 10 societies – Symbolic-operational ideological alignment ===\n")
top10_symbolic <- polar %>%
  arrange(desc(symbolic_op_align)) %>%
  slice(1:10) %>%
  dplyr::select(country, symbolic_op_align)
print(top10_symbolic)

cat("\n=== Table 2: Top 10 societies – Within-operational ideological alignment ===\n")
top10_within <- polar %>%
  arrange(desc(within_op_align)) %>%
  slice(1:10) %>%
  dplyr::select(country, within_op_align)
print(top10_within)


################################################################################
# SECTION 7 — MAIN TEXT FIGURE 4: Scatter Plots Between Manifestations
#
# Three panels (A–C) show pairwise scatter plots among the three manifestations,
# with Spearman ρ annotations.
################################################################################

# Panel A: Disagreement vs. Symbolic-operational ideological alignment
f4a <- ggplot(polar, aes(x = disagreement, y = symbolic_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Disagreement",
       y = "Symbolic-operational ideological alignment") +
  annotate("text", x = Inf, y = Inf,
           label = "Spearman \u03c1 = \u22120.387 (p < .001)",
           hjust = 1.1, vjust = 1.5, size = 4)
f4a <- ggdraw(f4a) +
  draw_plot_label("A", x = 0, y = 1, hjust = -0.3, vjust = 1.5,
                  size = 20, fontface = "bold")

# Panel B: Disagreement vs. Within-operational ideological alignment
f4b <- ggplot(polar, aes(x = disagreement, y = within_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Disagreement",
       y = "Within-operational ideological alignment") +
  annotate("text", x = Inf, y = Inf,
           label = "Spearman \u03c1 = \u22120.335 (p = .003)",
           hjust = 1.1, vjust = 1.5, size = 4)
f4b <- ggdraw(f4b) +
  draw_plot_label("B", x = 0, y = 1, hjust = -0.3, vjust = 1.5,
                  size = 20, fontface = "bold")

# Panel C: Symbolic-operational vs. Within-operational ideological alignment
f4c <- ggplot(polar, aes(x = symbolic_op_align, y = within_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Symbolic-operational ideological alignment",
       y = "Within-operational ideological alignment") +
  annotate("text", x = Inf, y = Inf,
           label = "Spearman \u03c1 = 0.571 (p < .001)",
           hjust = 1.1, vjust = 1.5, size = 4)
f4c <- ggdraw(f4c) +
  draw_plot_label("C", x = 0, y = 1, hjust = -0.3, vjust = 1.5,
                  size = 20, fontface = "bold")

Figure4 <- ggpubr::ggarrange(f4a, f4b, f4c, ncol = 1)
ggsave("Figure4_scatter_manifestations.png",
       plot = Figure4, height = 12, width = 8)
cat("Saved Figure4_scatter_manifestations.png\n")


################################################################################
# SECTION 8 — MAIN TEXT FIGURE 3: World Maps
#
# Three-panel choropleth map showing the global distribution of Disagreement
# (Panel A), Symbolic-operational ideological alignment (Panel B), and
# Within-operational ideological alignment (Panel C).
# Grey = NA (country not in sample).
################################################################################

# Load the Gleditsch–Ward world map (2015 borders)
world_mp <- cshp(date = as.Date("2015-01-01"), useGW = TRUE)
map2 <- world_mp %>%
  mutate(country_name = as.character(country_name))

# Harmonise country names in the shapefile to match those in polar$country
map2$country_name[map2$country_name == "Korea, Republic of"]               <- "South_Korea"
map2$country_name[map2$country_name == "Russia (Soviet Union)"]             <- "Russia"
map2$country_name[map2$country_name == "Belarus (Byelorussia)"]             <- "Belarus"
map2$country_name[map2$country_name == "Bosnia-Herzegovina"]                <- "Bosnia_Herzegovina"
map2$country_name[map2$country_name == "Czech Republic"]                    <- "Czechia"
map2$country_name[map2$country_name == "German Federal Republic"]           <- "Germany"
map2$country_name[map2$country_name == "Italy/Sardinia"]                    <- "Italy"
map2$country_name[map2$country_name == "Cote D'Ivoire"]                     <- "Ivory Coast"
map2$country_name[map2$country_name == "Kyrgyz Republic"]                   <- "Kyrgyzstan"
map2$country_name[map2$country_name == "Macedonia (FYROM/North Macedonia)"] <- "North_Macedonia"
map2$country_name[map2$country_name == "Zimbabwe (Rhodesia)"]               <- "Zimbabwe"
map2$country_name[map2$country_name == "Rumania"]                           <- "Romania"
map2$country_name[map2$country_name == "United States of America"]          <- "United_States"
map2$country_name[map2$country_name == "United Kingdom"]                    <- "Great_Britain"
map2$country_name[map2$country_name == "New Zealand"]                       <- "New_Zealand"

# Exact-match merge (max_dist = 0 disables fuzzy matching)
map3 <- stringdist_left_join(map2, polar,
                              by = c("country_name" = "country"),
                              method = "jw", max_dist = 0) %>%
  sf::st_as_sf()

# Panel A: Disagreement
f3a <- ggplot(map3) +
  geom_sf(aes(fill = disagreement), color = NA) +
  scale_fill_gradientn(
    colours  = rev(viridis::magma(4)),
    name     = "",
    na.value = "grey90",
    breaks   = c(0.3, 0.4, 0.5, 0.6),
    limits   = c(0.3, 0.6),
    oob      = scales::squish) +
  theme_void() +
  labs(title = "A. Disagreement") +
  theme(plot.title = element_text(hjust = 0.5, size = 12),
        legend.position = "right")

# Panel B: Symbolic-operational ideological alignment
f3b <- ggplot(map3) +
  geom_sf(aes(fill = symbolic_op_align), color = NA) +
  scale_fill_gradientn(
    colours  = rev(viridis::magma(5)),
    name     = "",
    na.value = "grey90",
    breaks   = c(0, 0.5, 1, 1.5),
    limits   = c(0, 1.5),
    oob      = scales::squish) +
  theme_void() +
  labs(title = "B. Symbolic-operational ideological alignment") +
  theme(plot.title = element_text(hjust = 0.5, size = 12),
        legend.position = "right")

# Panel C: Within-operational ideological alignment
f3c <- ggplot(map3) +
  geom_sf(aes(fill = within_op_align), color = NA) +
  scale_fill_gradientn(
    colours  = rev(viridis::magma(5)),
    name     = "",
    na.value = "grey90",
    breaks   = c(5, 7.5, 10, 12.5, 15),
    limits   = c(4, 15),
    oob      = scales::squish) +
  theme_void() +
  labs(title = "C. Within-operational ideological alignment") +
  theme(plot.title = element_text(hjust = 0.5, size = 12),
        legend.position = "right")

Figure3 <- ggpubr::ggarrange(f3a, f3b, f3c, ncol = 1)
ggsave("Figure3_maps.png", plot = Figure3, height = 12, width = 10)
cat("Saved Figure3_maps.png\n")


################################################################################
# SECTION 9 — PREPARE COUNTRY-LEVEL COVARIATES FOR REGRESSIONS
#
# Merge the polarization scores with country-level predictor variables.
# Predictors: HDI, GDP per capita (log), political stability, political
# parallelism, total population (log), press freedom, ethnic fractionalization.
# Source: countrycode.xlsx Sheet 2, supplemented with V-Dem v14.
################################################################################

# 9a. Load V-Dem data ──────────────────────────────────────────────────────────
vdem <- readRDS("V-Dem-CY-Full+Others-v14.rds")

vdem1 <- vdem %>%
  filter(year >= 2017) %>%
  dplyr::select(
    country_name, country_text_id, year,
    v2x_polyarchy,    # Electoral democracy index (not used in main regressions)
    v2mebias,         # Media bias toward gov. party = Political parallelism
    v2xps_party,      # Party System Institutionalization (PSI; SI Sec E)
    v2pscomprg,       # Party competition across regions (SI Sec E)
    v2psplats,        # Distinct party platforms (SI Sec E)
    v2smpolsoc        # Issue polarization (reversed-coded; SI Sec E)
  ) %>%
  dplyr::rename(
    polparrel     = v2mebias,      # Political parallelism
    PSI           = v2xps_party,   # Party system institutionalization
    party_compr   = v2pscomprg,    # Party competition across regions
    party_platf   = v2psplats,     # Distinct party platforms
    issuepolar    = v2smpolsoc     # Issue polarization (V-Dem; SI Sec E)
  )

# 9b. Load country-level covariates from countrycode.xlsx Sheet 2 ─────────────
countrydat_raw <- read_excel("countrycode.xlsx", sheet = "Sheet2")

# Merge with V-Dem variables and average across survey years within each country
countrydat <- merge(countrydat_raw, vdem1,
                    by = c("country_text_id", "year"), all.x = TRUE) %>%
  group_by(country) %>%
  summarise(across(c(polstability, population, HDI, pressfreedom,
                     ethnicfract, GDP, polparrel),
                   ~ mean(.x, na.rm = TRUE))) %>%
  ungroup()

# 9c. Merge with polarization scores ──────────────────────────────────────────
countrydat <- merge(countrydat, polar, all.x = TRUE, by = "country")

# 9d. Log-transform skewed variables ──────────────────────────────────────────
countrydat$population <- log(countrydat$population)
countrydat$GDP        <- log(countrydat$GDP)

# 9e. Listwise deletion and z-score standardisation ───────────────────────────
# The final regression sample is 67 societies (Table C2 in SI).
countrydat_analysis  <- countrydat %>% na.omit()
countrydat_scaled    <- countrydat_analysis %>%
  mutate(across(where(is.numeric), scale)) %>%
  as.data.frame()
# Add squared terms for potential curvilinear effects (see Sec 10)
countrydat_scaled$HDI2 <- (countrydat_scaled$HDI)^2
countrydat_scaled$GDP2 <- (countrydat_scaled$GDP)^2

cat("\nRegression sample: N =", nrow(countrydat_scaled), "societies\n")


################################################################################
# SECTION 10 — MAIN TEXT FIGURE 5: OLS Regression Coefficient Plot
#
# Three OLS models predicting each manifestation from HDI (linear + squared),
# political stability, total population (log), press freedom, ethnic
# fractionalization, and political parallelism. Coefficients are standardised.
# Results are visualised as a faceted dot-whisker plot (95% CIs).
################################################################################

# 10a. Fit the three OLS models ───────────────────────────────────────────────
M_disagree  <- lm(disagreement      ~ HDI + HDI2 + polstability + population +
                    pressfreedom + ethnicfract + polparrel,
                  data = countrydat_scaled)

M_symbolic  <- lm(symbolic_op_align ~ HDI + HDI2 + polstability + population +
                    pressfreedom + ethnicfract + polparrel,
                  data = countrydat_scaled)

M_within    <- lm(within_op_align   ~ HDI + HDI2 + polstability + population +
                    pressfreedom + ethnicfract + polparrel,
                  data = countrydat_scaled)

# Print regression table
screenreg(
  list(M_disagree, M_symbolic, M_within),
  digits = 3,
  stars  = c(0.1, 0.05, 0.01, 0.001),
  custom.model.names = c("Disagreement",
                         "Symbolic-operational alignment",
                         "Within-operational alignment")
)

# 10b. Build the coefficient plot data frame ──────────────────────────────────
extract_coefs <- function(model, dv_label) {
  coefs <- summary(model)$coefficients
  data.frame(
    coef       = coefs[, 1],
    se         = coefs[, 2],
    p_value    = coefs[, 4],
    predictor  = rownames(coefs),
    dv         = dv_label,
    stringsAsFactors = FALSE
  )
}

plotdf <- rbind(
  extract_coefs(M_disagree, "Disagreement"),
  extract_coefs(M_symbolic, "Symbolic-operational\nideological alignment"),
  extract_coefs(M_within,   "Within-operational\nideological alignment")
) %>%
  filter(predictor != "(Intercept)")

# Significance stars
plotdf$sig <- ifelse(plotdf$p_value < 0.001, "***",
              ifelse(plotdf$p_value < 0.01,  "**",
              ifelse(plotdf$p_value < 0.05,  "*",
              ifelse(plotdf$p_value < 0.1,   "^", ""))))
plotdf$label <- paste0(round(plotdf$coef, 3), plotdf$sig)

# Rename and reorder predictors for display
pred_labels <- c(
  "HDI"          = "HDI",
  "HDI2"         = "HDI (squared)",
  "polstability" = "Political Stability",
  "polparrel"    = "Political Parallelism",
  "population"   = "Total Population",
  "pressfreedom" = "Press Freedom",
  "ethnicfract"  = "Ethnic Fractionalization"
)
plotdf$predictor <- pred_labels[plotdf$predictor]
pred_order <- rev(c("HDI", "HDI (squared)", 
                    "Political Stability",
                    "Political Parallelism",
                    "Total Population", 
                    "Press Freedom",
                    "Ethnic Fractionalization"))
plotdf$predictor <- factor(plotdf$predictor, levels = pred_order)
plotdf$dv <- factor(plotdf$dv, levels = c("Disagreement",
                                           "Symbolic-operational\nideological alignment",
                                           "Within-operational\nideological alignment"))

# 10c. Plot ───────────────────────────────────────────────────────────────────
Figure5 <- ggplot(plotdf, aes(x = predictor, y = coef, color = dv)) +
  facet_wrap(~dv, scales = "free_x") +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, hjust = 0.5,
            size = 3, color = "black") +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "") +
  theme_bw() +
  theme(
    axis.text.y  = element_text(size = 10),
    axis.text.x  = element_text(size = 8),
    legend.position = "none",
    strip.background = element_rect(fill = "grey90"),
    strip.text = element_text(size = 9)
  )
ggsave("Figure5_coefplot.png", plot = Figure5, height = 7, width = 11)
cat("Saved Figure5_coefplot.png\n")


################################################################################
# SECTION 11 — MAIN TEXT FIGURE 6: Scatter Plots — HDI vs. Each Manifestation
#
# Three panels (A–C) showing the bivariate (quadratic) relationship between
# HDI and each polarization manifestation.
################################################################################

f6a <- ggplot(countrydat_analysis, aes(x = HDI, y = disagreement)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Disagreement", x = "HDI")
f6a <- ggdraw(f6a) +
  draw_plot_label("A", x = 0, y = 0.97, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

f6b <- ggplot(countrydat_analysis, aes(x = HDI, y = symbolic_op_align)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Symbolic-operational ideological alignment", x = "HDI")
f6b <- ggdraw(f6b) +
  draw_plot_label("B", x = 0, y = 1, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

f6c <- ggplot(countrydat_analysis, aes(x = HDI, y = within_op_align)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Within-operational ideological alignment", x = "HDI")
f6c <- ggdraw(f6c) +
  draw_plot_label("C", x = 0, y = 1, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

Figure6 <- ggpubr::ggarrange(f6a, f6b, f6c, ncol = 1)
ggsave("Figure6_HDI_scatter.png", plot = Figure6, height = 12, width = 8)
cat("Saved Figure6_HDI_scatter.png\n")


################################################################################
# SECTION 12 — MAIN TEXT TABLE 3: OLS Models Predicting Disagreement
#
# Four models: (1) symbolic-op alignment only; (2) symbolic-op + HDI controls;
# (3) within-op alignment only; (4) within-op + HDI controls.
# Tests whether the negative inter-manifestation correlations are explained
# by modernization (HDI) rather than by a direct relationship.
################################################################################

T3_m1 <- lm(disagreement ~ symbolic_op_align,
             data = countrydat_scaled)
T3_m2 <- lm(disagreement ~ symbolic_op_align + HDI + HDI2,
             data = countrydat_scaled)
T3_m3 <- lm(disagreement ~ within_op_align,
             data = countrydat_scaled)
T3_m4 <- lm(disagreement ~ within_op_align + HDI + HDI2,
             data = countrydat_scaled)

cat("\n=== Table 3: OLS Models Predicting Disagreement ===\n")
screenreg(list(T3_m1, T3_m2, T3_m3, T3_m4), digits = 3)


################################################################################
# SECTION 13 — MAIN TEXT FIGURE 1 / SI FIGURE H1: Belief Network Simulation
#
# Illustrates how naïve pairwise correlations produce spurious edges, whereas
# the LASSO-regularised partial-correlation network (EBICglasso) recovers the
# true chain structure.
# This code produces both Figure 1 in the main text and Figure H1 in the SI.
################################################################################

library(MASS)
set.seed(123)

n_sim <- 1000   # simulated individuals
p_sim <- 8      # attitudinal variables

# Specify a precision matrix (inverse covariance) with a chain structure:
# only adjacent pairs (A–B, B–C, …, G–H) have direct conditional dependencies.
precision <- diag(p_sim)                       # identity matrix
for (i in 1:(p_sim - 1)) {
  precision[i, i + 1] <- precision[i + 1, i] <- 0.3
}

# Generate data from the multivariate normal implied by this precision matrix
cov_mat  <- solve(precision)
sim_data <- mvrnorm(n = n_sim, mu = rep(0, p_sim), Sigma = cov_mat)
colnames(sim_data) <- LETTERS[1:p_sim]

# Panel A: Pairwise Pearson correlation network
# Many non-adjacent edges appear due to transitive (indirect) associations.
png("Figure1_H1_panel_A_correlation.png", width = 600, height = 600)
qgraph(cor(sim_data),
       layout     = "circle",
       color      = "white",
       labels     = LETTERS[1:p_sim],
       edge.color = "black",
       title      = "Panel A: Correlation Network")
dev.off()

# Panel B: Regularised partial-correlation network (EBICglasso)
# Recovers the true chain structure by conditioning on all other variables.
pcor_mat <- qgraph::cor_auto(sim_data)
pcor_net <- qgraph::EBICglasso(pcor_mat, n = n_sim)
png("Figure1_H1_panel_B_partial_correlation.png", width = 600, height = 600)
qgraph(pcor_net,
       layout     = "circle",
       color      = "white",
       labels     = LETTERS[1:p_sim],
       edge.color = "black",
       title      = "Panel B: Partial Correlation Network (EBICglasso)")
dev.off()

cat("Saved Figure1_H1 panels.\n")


################################################################################
# ─────────────────────────────────────────────────────────────────────────────
#  SUPPLEMENTARY INFORMATION
# ─────────────────────────────────────────────────────────────────────────────
################################################################################


################################################################################
# SI SECTION A — Table A1: Full Item List with Disagree–Overall Correlations

# Load item-level disagreement matrix (computed in Section 4 and saved to CSV)
# OR re-compute it here:
# item_disagree_df <- as.data.frame(sapply(countrylist, disagreement_itemlevel_fun))
# For each item: correlation between its column of country-level scores and
# the vector of overall country-level disagreement scores.

item_disagree_mat <- read.csv("item_disagreement.csv",
                               row.names = "item",
                               stringsAsFactors = FALSE)
# item_disagree_mat: rows = items (34), columns = countries

# Overall disagreement scores (matching country order)
overall_disagree <- polar$disagreement[match(colnames(item_disagree_mat),
                                              polar$country)]

# Spearman ρ for each item
table_A1_rho <- sapply(rownames(item_disagree_mat), function(item) {
  item_scores <- as.numeric(item_disagree_mat[item, ])
  cor(item_scores, overall_disagree, method = "spearman",
      use = "pairwise.complete.obs")
})

cat("\n=== SI Table A1: Spearman ρ (item disagree ~ overall disagree) ===\n")
print(round(sort(table_A1_rho, decreasing = TRUE), 3))


################################################################################
# SI SECTION B — Excluding Additional Ideological Items (Figure B1)
#
# The "governmental/democratic" items (E114–E117, E224–E229, E233, E235) may
# conflate ideological preferences with orientations toward democracy per se.
# We re-run the main analyses excluding these 11 items and keep only the 23
# economic and social items.
#
# NOTE: This section re-runs the full network estimation on the reduced item set.
#       Estimated additional run time: 15–45 minutes.
################################################################################

# Governmental/democratic items to exclude
govt_items <- c("strongleader", "experts", "dompol",
                 "taxrich", "relig", "election", "aidunempl",
                 "army", "civilright", "womenright", "impdem")

issuelist_B <- setdiff(issuelist, govt_items)   # 23 remaining items
cat("\nSection B: reduced item set has", length(issuelist_B), "items\n")

# Re-run estimation on reduced item set ────────────────────────────────────────
# (temporarily redefine helper functions to use issuelist_B)
within_op_fun_B    <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist_B, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    g   <- estimateNetwork(dat, default = "mgm")
    adj <- abs(g$graph)
    sum(adj[upper.tri(adj)])
  }, error = function(e) NA_real_)
}

symbolic_op_fun_B  <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, c("ideology", issuelist_B), drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    g   <- estimateNetwork(dat, default = "mgm")
    ct  <- centrality(g)
    ct[["InDegree"]][1]
  }, error = function(e) NA_real_)
}

disagreement_fun_B <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, issuelist_B, drop = FALSE]
    dat <- labelled::remove_labels(dat)
    item_disagree <- function(item_name) {
      x <- dat[[item_name]]
      x <- x[!is.na(x)]
      if (length(x) < 10) return(NA_real_)
      1 - agrmt::consensus(as.vector(round(
        questionr::wtd.table(x, weights = rep(1, length(x))), 0)))
    }
    mean(sapply(issuelist_B, item_disagree), na.rm = TRUE)
  }, error = function(e) NA_real_)
}

cat("Section B: running network estimation on reduced item set...\n")
polar_B <- data.frame(
  country           = countrylist,
  disagreement      = sapply(countrylist, disagreement_fun_B),
  symbolic_op_align = sapply(countrylist, symbolic_op_fun_B),
  within_op_align   = sapply(countrylist, within_op_fun_B),
  stringsAsFactors  = FALSE
) %>% filter(!is.na(disagreement) & !is.na(symbolic_op_align) & !is.na(within_op_align))
# write.csv(polar_B, "dimension_summary_B_reduced_items.csv", row.names = FALSE)

# Inter-manifestation correlations with reduced item set
cat("\nSection B: inter-manifestation correlations (reduced items)\n")
cat("Disagreement ~ Symbolic-op:\n")
print(cor.test(polar_B$disagreement, polar_B$symbolic_op_align, method = "spearman"))
cat("Disagreement ~ Within-op:\n")
print(cor.test(polar_B$disagreement, polar_B$within_op_align,   method = "spearman"))
cat("Symbolic-op ~ Within-op:\n")
print(cor.test(polar_B$symbolic_op_align, polar_B$within_op_align, method = "spearman"))

# Figure B1: regression coefficient plot with reduced items ───────────────────
countrydat_B <- merge(countrydat_analysis[, !names(countrydat_analysis) %in%
                        c("disagreement","symbolic_op_align","within_op_align")],
                      polar_B, all.x = TRUE, by = "country") %>% na.omit()
countrydat_B_scaled <- countrydat_B %>%
  mutate(across(where(is.numeric), scale)) %>%
  as.data.frame()
countrydat_B_scaled$HDI2 <- (countrydat_B_scaled$HDI)^2

MB1 <- lm(disagreement      ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_B_scaled)
MB2 <- lm(symbolic_op_align ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_B_scaled)
MB3 <- lm(within_op_align   ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_B_scaled)

plotdf_B <- rbind(
  extract_coefs(MB1, "Disagreement"),
  extract_coefs(MB2, "Symbolic-operational\nideological alignment"),
  extract_coefs(MB3, "Within-operational\nideological alignment")
) %>% filter(predictor != "(Intercept)")

plotdf_B$sig <- ifelse(plotdf_B$p_value < 0.001, "***",
               ifelse(plotdf_B$p_value < 0.01,  "**",
               ifelse(plotdf_B$p_value < 0.05,  "*",
               ifelse(plotdf_B$p_value < 0.1,   "^", ""))))
plotdf_B$label    <- paste0(round(plotdf_B$coef, 3), plotdf_B$sig)
plotdf_B$predictor <- pred_labels[plotdf_B$predictor]
plotdf_B$predictor <- factor(plotdf_B$predictor, levels = pred_order)
plotdf_B$dv <- factor(plotdf_B$dv, levels = c("Disagreement",
                                                "Symbolic-operational\nideological alignment",
                                                "Within-operational\nideological alignment"))

FigureB1 <- ggplot(plotdf_B, aes(x = predictor, y = coef, color = dv)) +
  facet_wrap(~dv, scales = "free_x") +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, hjust = 0.5,
            size = 3, color = "black") +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "",
       subtitle = "SI Section B: Reduced item set (economic + social items only, N = 23)") +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10), legend.position = "none",
        strip.background = element_rect(fill = "grey90"),
        strip.text = element_text(size = 9))
ggsave("SI_Figure_B1_coefplot_reduced_items.png",
       plot = FigureB1, height = 5, width = 9)
cat("Saved SI_Figure_B1_coefplot_reduced_items.png\n")


################################################################################
# SI SECTION C — Country-Level Variable Descriptives (Tables C1–C2)
#
# Table C1 is a narrative description table (see SI document).
# Table C2 provides descriptive statistics before z-score standardisation.
################################################################################

cat("\n=== SI Table C2: Descriptive statistics (before standardisation) ===\n")
desc_vars <- c("HDI", "GDP", "polstability", "polparrel",
               "population", "pressfreedom", "ethnicfract",
               "symbolic_op_align", "within_op_align", "disagreement")
print(psych::describe(countrydat_analysis[, desc_vars], fast = TRUE))


################################################################################
# SI SECTION D — Country-Year as the Unit of Analysis (Figure D1)
#
# To address temporal variation (some countries surveyed in multiple years),
# we treat each country-year combination as a separate unit. We re-run the
# network estimation by country-year, yielding N = 96 observations for
# the regression analyses.
#
# NOTE: This requires a second full pass of network estimation. Estimated
#       run time: 60–120 minutes (more units, smaller per-unit samples).
################################################################################

# 9b re-do: use country-year combinations rather than aggregating ──────────────
countrydat_cy <- merge(countrydat_raw, vdem1,
                        by = c("country_text_id", "year"), all.x = TRUE)
countrydat_cy <- countrydat_cy[,c(2, 3, 10,11, 13, 14, 15, 18)]

# Merge with country-year covariates
countrydat_cy2 <- merge(countrydat_cy, polar,
                         by = "country", all.x = TRUE) %>%
  na.omit()
countrydat_cy2$population <- log(countrydat_cy2$population)
cat("Country-year regression N =", nrow(countrydat_cy2), "\n")

countrydat_cy_scaled <- countrydat_cy2 %>%
  mutate(across(where(is.numeric), scale)) %>%
  as.data.frame()
countrydat_cy_scaled$HDI2 <- (countrydat_cy_scaled$HDI)^2

MD1 <- lm(disagreement      ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_cy_scaled)
MD2 <- lm(symbolic_op_align ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_cy_scaled)
MD3 <- lm(within_op_align   ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_cy_scaled)

plotdf_D <- rbind(
  extract_coefs(MD1, "Disagreement"),
  extract_coefs(MD2, "Symbolic-operational\nideological alignment"),
  extract_coefs(MD3, "Within-operational\nideological alignment")
) %>% filter(predictor != "(Intercept)")

plotdf_D$sig <- ifelse(plotdf_D$p_value < 0.001, "***",
               ifelse(plotdf_D$p_value < 0.01,  "**",
               ifelse(plotdf_D$p_value < 0.05,  "*",
               ifelse(plotdf_D$p_value < 0.1,   "^", ""))))
plotdf_D$label    <- paste0(round(plotdf_D$coef, 3), plotdf_D$sig)
plotdf_D$predictor <- pred_labels[plotdf_D$predictor]
plotdf_D$predictor <- factor(plotdf_D$predictor, levels = pred_order)
plotdf_D$dv <- factor(plotdf_D$dv, levels = c("Disagreement",
                                                "Symbolic-operational\nideological alignment",
                                                "Within-operational\nideological alignment"))

FigureD1 <- ggplot(plotdf_D, aes(x = predictor, y = coef, color = dv)) +
  facet_wrap(~dv, scales = "free_x") +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, hjust = 0.5,
            size = 3, color = "black") +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "",
       subtitle = paste0("SI Section D: Country-year as unit of analysis (N = ",
                          nrow(countrydat_cy2), ")")) +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10), legend.position = "none",
        strip.background = element_rect(fill = "grey90"),
        strip.text = element_text(size = 9))
ggsave("SI_Figure_D1_coefplot_country_year.png",
       plot = FigureD1, height = 5, width = 9)
cat("Saved SI_Figure_D1_coefplot_country_year.png\n")


################################################################################
# SI SECTION E — Measurement Validation (Figures E1–E2)
#
# Figure E1: Correlation matrix among all eight alternative measures
#   (three for disagreement, two for symbolic-op, three for within-op).
#   Shows construct validity: same-construct measures correlate positively;
#   measures of different constructs correlate negatively.
#
# Figure E2: Correlations between each manifestation and five theoretically
#   related V-Dem/WGI indicators (political stability, issue polarization,
#   party system institutionalization, party competition, distinct platforms).
################################################################################

# Figure E1 ───────────────────────────────────────────────────────────────────
# Column names must match those in dimension_summary.csv
measures_E1 <- polar %>%
  dplyr::select(
    `Disagreement\n(ordinal dispersion)` = disagreement,
    `Disagreement\n(bimodality coef.)`   = bimodality_coef,
    `Disagreement\n(mean variance)`      = avg_variance,
    `Symbolic-op\n(InDegree)`            = symbolic_op_align,
    `Symbolic-op\n(pairwise cor.)`       = avg_cor_ideology,
    `Within-op\n(sum of weights)`        = within_op_align,
    `Within-op\n(network density)`       = density,
    `Within-op\n(pairwise cor.)`         = avg_pairwise_cor
  ) %>%
  as.matrix()

cor_E1 <- cor(measures_E1, method = "spearman", use = "pairwise.complete.obs")

png("SI_Figure_E1_measure_correlations.png", width = 500, height = 500)

corrplot(
  cor_E1,
  method  = "number",
  tl.col  = "black",
  tl.srt  = 90,
  tl.cex  = 01)
dev.off()
cat("Saved SI_Figure_E1_measure_correlations.png\n")

# Figure E2 ───────────────────────────────────────────────────────────────────
# Merge manifestation scores with V-Dem validation indicators.
# V-Dem variables needed: PSI (v2xps_party), party_compr (v2pscomprg),
# party_platf (v2psplats), issuepolar (v2smpolsoc); polstability.
countrydat <- read_excel("countrycode.xlsx", sheet = "Sheet2")
countrydat <- merge(countrydat, vdem1, by=c("country_text_id", "year"), all.x=T)
countrydat <- countrydat %>%
  group_by(country)%>%
  summarise(across(c(polstability, population,
                     HDI, pressfreedom, ethnicfract, 
                     polparrel, GDP, PSI, party_compr, issuepolar, party_platf), 
                   ~ mean(.x, na.rm = TRUE)))

countrydat <- merge(countrydat, polar, all.x=T, by ="country")

validation_dat <- merge(countrydat,
                        polar[, c("country","disagreement",
                                   "symbolic_op_align","within_op_align")],
                         by = "country", all.x = TRUE) %>%
  na.omit()

cor_E2 <- cor(validation_dat[,c(2,11,9,10,12:15)], method = "spearman",
              use = "pairwise.complete.obs")

png("SI_Figure_E2_external_validation.png", width = 500, height = 300)
corrplot::corrplot(
  cor_E2[6:8, 1:5],
  method  = "number",
  tl.col  = "black",
  tl.cex  = 0.9)
dev.off()
cat("Saved SI_Figure_E2_external_validation.png\n")


################################################################################
# SI SECTION F — Within-operational Alignment Measured by Network Density
#
# Replicates Figure 4 (inter-manifestation scatter), Figure 3 (world map),
# and Figure 5 (regression coefficient plot) using network density instead
# of sum of edge weights as the within-operational ideological alignment measure.
################################################################################

# Use the density column already computed in polar
polar_F <- polar %>%
  rename(within_op_align_density = density)

# Figure F1: Correlation matrix using density
cat("\n=== SI Section F: inter-manifestation correlations (density measure) ===\n")
cat("Disagreement ~ density:\n")
print(cor.test(polar_F$disagreement, polar_F$within_op_align_density,
               method = "spearman"))
cat("Symbolic-op ~ density:\n")
print(cor.test(polar_F$symbolic_op_align, polar_F$within_op_align_density,
               method = "spearman"))


# Figure F2: World map using network density
map3_F <- stringdist_left_join(map2, polar_F,
                                by = c("country_name" = "country"),
                                method = "jw", max_dist = 0) %>%
  sf::st_as_sf()
FigureF2 <- ggplot(map3_F) +
  geom_sf(aes(fill = within_op_align_density), color = NA) +
  scale_fill_gradientn(colours = rev(viridis::magma(5)), name = "",
                        na.value = "grey90", oob = scales::squish) +
  theme_void() +
  labs(title = "SI Figure F2: Within-operational alignment (network density)") +
  theme(plot.title = element_text(hjust = 0.5, size = 12))
ggsave("SI_Figure_F2_density_map.png", plot = FigureF2, height = 5, width = 10)

# Figure F3: Regression with density as DV
countrydat_F <- merge(countrydat_scaled,
                       polar_F[, c("country","within_op_align_density")],
                       by = "country", all.x = TRUE) %>% na.omit()
countrydat_F_scaled <- countrydat_F %>%
  mutate(across(where(is.numeric), scale)) %>%
  as.data.frame()
countrydat_F_scaled$HDI2 <- (countrydat_F_scaled$HDI)^2

MF3 <- lm(within_op_align_density ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_F_scaled)
cat("\n=== SI Section F: regression with network density as DV ===\n")
texreg::screenreg(list(MF3), digits = 3,
                  custom.model.names = "Within-op alignment (density)")

plotdf_F3 <- extract_coefs(MF3, "Within-op alignment (density)") %>%
  filter(predictor != "(Intercept)")
plotdf_F3$sig      <- ifelse(plotdf_F3$p_value < 0.05, "*", "")
plotdf_F3$label    <- paste0(round(plotdf_F3$coef, 3), plotdf_F3$sig)
plotdf_F3$predictor <- pred_labels[plotdf_F3$predictor]
plotdf_F3$predictor <- factor(plotdf_F3$predictor, levels = pred_order)
FigureF3 <- ggplot(plotdf_F3, aes(x = predictor, y = coef)) +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, size = 3) +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "",
       subtitle = "SI Figure F3: DV = Within-op alignment (network density)") +
  theme_bw() + theme(axis.text.y = element_text(size = 10))
ggsave("SI_Figure_F3_density_coefplot.png", plot = FigureF3, height = 5, width = 7)
cat("Saved SI Section F figures.\n")


################################################################################
# SI SECTION G — Modernization Measured by GDP per Capita (Figures G1–G2)
################################################################################

MG1 <- lm(disagreement      ~ GDP + GDP2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_scaled)
MG2 <- lm(symbolic_op_align ~ GDP + GDP2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_scaled)
MG3 <- lm(within_op_align   ~ GDP + GDP2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_scaled)

cat("\n=== SI Section G: regressions with GDP per capita ===\n")
texreg::screenreg(list(MG1, MG2, MG3), digits = 3,
                  custom.model.names = c("Disagreement",
                                         "Symbolic-op alignment",
                                         "Within-op alignment"))

# Figure G1: coefficient plot (GDP version)
pred_labels_G <- c(
  "GDP"          = "GDP per capita (log)",
  "GDP2"         = "GDP per capita (log, squared)",
  "polstability" = "Political Stability",
  "population"   = "Total Population (log)",
  "pressfreedom" = "Press Freedom",
  "ethnicfract"  = "Ethnic Fractionalization",
  "polparrel"    = "Political Parallelism"
)
pred_order_G <- rev(c("GDP per capita (log)", "GDP per capita (log, squared)",
                       "Political Stability", "Total Population (log)",
                       "Press Freedom", "Ethnic Fractionalization",
                       "Political Parallelism"))

plotdf_G <- rbind(
  extract_coefs(MG1, "Disagreement"),
  extract_coefs(MG2, "Symbolic-operational\nideological alignment"),
  extract_coefs(MG3, "Within-operational\nideological alignment")
) %>% filter(predictor != "(Intercept)")

plotdf_G$sig <- ifelse(plotdf_G$p_value < 0.001, "***",
               ifelse(plotdf_G$p_value < 0.01,  "**",
               ifelse(plotdf_G$p_value < 0.05,  "*",
               ifelse(plotdf_G$p_value < 0.1,   "^", ""))))
plotdf_G$label    <- paste0(round(plotdf_G$coef, 3), plotdf_G$sig)
plotdf_G$predictor <- pred_labels_G[plotdf_G$predictor]
plotdf_G$predictor <- factor(plotdf_G$predictor, levels = pred_order_G)
plotdf_G$dv <- factor(plotdf_G$dv, levels = c("Disagreement",
                                                "Symbolic-operational\nideological alignment",
                                                "Within-operational\nideological alignment"))

FigureG1 <- ggplot(plotdf_G, aes(x = predictor, y = coef, color = dv)) +
  facet_wrap(~dv, scales = "free_x") +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, hjust = 0.5,
            size = 3, color = "black") +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "",
       subtitle = "SI Figure G1: GDP per capita as modernization measure") +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10), legend.position = "none",
        strip.background = element_rect(fill = "grey90"),
        strip.text = element_text(size = 9))
ggsave("SI_Figure_G1_GDP_coefplot.png", plot = FigureG1, height = 5, width = 9)

# Figure G2: Scatter plots — GDP vs. each manifestation
fG2a <- ggplot(countrydat_analysis, aes(x = GDP, y = disagreement)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Disagreement", x = "GDP per capita (log-transformed)")
fG2a <- ggdraw(fG2a) +
  draw_plot_label("A", x = 0, y = 1, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

fG2b <- ggplot(countrydat_analysis, aes(x = GDP, y = symbolic_op_align)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Symbolic-operational ideological alignment",
       x = "GDP per capita (log-transformed)")
fG2b <- ggdraw(fG2b) +
  draw_plot_label("B", x = 0, y = 1, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

fG2c <- ggplot(countrydat_analysis, aes(x = GDP, y = within_op_align)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2),
              colour = "black", linewidth = 0.8, se = FALSE) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 2.5) +
  theme_bw() +
  labs(y = "Within-operational ideological alignment",
       x = "GDP per capita (log-transformed)")
fG2c <- ggdraw(fG2c) +
  draw_plot_label("C", x = 0, y = 1, hjust = -0.5, vjust = 1.5,
                  size = 20, fontface = "bold")

FigureG2 <- ggpubr::ggarrange(fG2a, fG2b, fG2c, ncol = 1)
ggsave("SI_Figure_G2_GDP_scatter.png", plot = FigureG2, height = 12, width = 8)
cat("Saved SI Section G figures.\n")


################################################################################
# SI SECTION H — Belief Network Simulation Detail (Figure H1)
#
# Same code as Section 13 (Main Figure 1), with additional descriptive text.
# Produces a side-by-side figure panel (Panel A: correlation network;
# Panel B: partial correlation network).
################################################################################

# See Section 13 above — Main Figure 1 / SI Figure H1 panels are produced
# there and saved as Figure1_H1_panel_A_correlation.png and
# Figure1_H1_panel_B_partial_correlation.png.
cat("SI Figure H1 panels were saved in Section 13.\n")
cat("File: Figure1_H1_panel_A_correlation.png\n")
cat("File: Figure1_H1_panel_B_partial_correlation.png\n")


################################################################################
# SI SECTION I — Weighted Results (Figures I1–I2, Table I1)
#
# Re-estimates all models using the EVS/WVS population weight (gwght) to
# correct for age, sex, education, and regional representativeness.
# The weighted disagreement measure uses wtd.table() with gwght as weights.
#
# NOTE: The network estimation (symbolic-op and within-op alignment) requires
#       the bootnet::estimateNetwork() 'weights' argument. This is computationally
#       intensive. Weighted results are reported for completeness; the unweighted
#       analyses are in the main text and confirmed robust here.
################################################################################

# Weighted disagreement function ──────────────────────────────────────────────
disagreement_wtd_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, c(issuelist, "gwght"), drop = FALSE]
    dat <- labelled::remove_labels(dat)

    item_disagree_wtd <- function(item_name) {
      x <- dat[[item_name]]
      w <- dat[["gwght"]]
      ok <- !is.na(x) & !is.na(w)
      x  <- x[ok]; w <- w[ok]
      if (length(x) < 10) return(NA_real_)
      freq_table <- as.vector(round(
        questionr::wtd.table(x, weights = w), 0))
      1 - agrmt::consensus(freq_table)
    }

    mean(sapply(issuelist, item_disagree_wtd), na.rm = TRUE)
  }, error = function(e) NA_real_)
}

# Weighted within-op and symbolic-op functions ─────────────────────────────────
within_op_wtd_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, c(issuelist, "gwght"), drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    w   <- dat$gwght
    dat <- dat[, issuelist]
    g   <- estimateNetwork(dat, default = "mgm", weights = w)
    adj <- abs(g$graph)
    sum(adj[upper.tri(adj)])
  }, error = function(e) NA_real_)
}

symbolic_op_wtd_fun <- function(n) {
  tryCatch({
    dat <- subset(wvs7, country == n)[, c("ideology", issuelist, "gwght"),
                                       drop = FALSE]
    dat <- labelled::remove_labels(dat)
    dat <- dat[complete.cases(dat), ]
    w   <- dat$gwght
    dat <- dat[, c("ideology", issuelist)]
    g   <- estimateNetwork(dat, default = "mgm", weights = w)
    ct  <- centrality(g)
    ct[["InDegree"]][1]
  }, error = function(e) NA_real_)
}

cat("Section I: running weighted network estimation...\n")
polar_I <- data.frame(
  country           = countrylist,
  disagreement      = sapply(countrylist, disagreement_wtd_fun),
  symbolic_op_align = sapply(countrylist, symbolic_op_wtd_fun),
  within_op_align   = sapply(countrylist, within_op_wtd_fun),
  stringsAsFactors  = FALSE
) %>%
  filter(!is.na(disagreement) & !is.na(symbolic_op_align) & !is.na(within_op_align))
write.csv(polar_I, "dimension_summary_I_weighted.csv", row.names = FALSE)

# Inter-manifestation correlations (weighted)
cat("\n=== SI Section I: inter-manifestation Spearman ρ (weighted) ===\n")
cor.test(polar_I$disagreement, polar_I$symbolic_op_align, method = "spearman")
cor.test(polar_I$disagreement, polar_I$within_op_align,   method = "spearman")
cor.test(polar_I$symbolic_op_align, polar_I$within_op_align, method = "spearman")

# Figure I1: scatter plots (weighted)
fI1a <- ggplot(polar_I, aes(x = disagreement, y = symbolic_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Disagreement (weighted)",
       y = "Symbolic-operational ideological alignment (weighted)")
fI1b <- ggplot(polar_I, aes(x = disagreement, y = within_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Disagreement (weighted)",
       y = "Within-operational ideological alignment (weighted)")
fI1c <- ggplot(polar_I, aes(x = symbolic_op_align, y = within_op_align)) +
  geom_point(size = 1.5) +
  ggrepel::geom_text_repel(aes(label = country), max.overlaps = 8, size = 3) +
  theme_bw() +
  labs(x = "Symbolic-operational alignment (weighted)",
       y = "Within-operational alignment (weighted)")
FigureI1 <- ggpubr::ggarrange(fI1a, fI1b, fI1c, ncol = 1)
ggsave("SI_Figure_I1_scatter_weighted.png",
       plot = FigureI1, height = 12, width = 8)

# Figure I2: regression coefficient plot (weighted)
countrydat_I <- merge(countrydat_analysis[, !names(countrydat_analysis) %in%
                        c("disagreement","symbolic_op_align","within_op_align")],
                       polar_I, by = "country", all.x = TRUE) %>% na.omit()
countrydat_I_scaled <- countrydat_I %>%
  mutate(across(where(is.numeric), scale)) %>%
  as.data.frame()
countrydat_I_scaled$HDI2 <- (countrydat_I_scaled$HDI)^2

MI1 <- lm(disagreement      ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_I_scaled)
MI2 <- lm(symbolic_op_align ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_I_scaled)
MI3 <- lm(within_op_align   ~ HDI + HDI2 + polstability + population +
             pressfreedom + ethnicfract + polparrel,
           data = countrydat_I_scaled)

plotdf_I <- rbind(
  extract_coefs(MI1, "Disagreement"),
  extract_coefs(MI2, "Symbolic-operational\nideological alignment"),
  extract_coefs(MI3, "Within-operational\nideological alignment")
) %>% filter(predictor != "(Intercept)")

plotdf_I$sig <- ifelse(plotdf_I$p_value < 0.001, "***",
               ifelse(plotdf_I$p_value < 0.01,  "**",
               ifelse(plotdf_I$p_value < 0.05,  "*",
               ifelse(plotdf_I$p_value < 0.1,   "^", ""))))
plotdf_I$label    <- paste0(round(plotdf_I$coef, 3), plotdf_I$sig)
plotdf_I$predictor <- pred_labels[plotdf_I$predictor]
plotdf_I$predictor <- factor(plotdf_I$predictor, levels = pred_order)
plotdf_I$dv <- factor(plotdf_I$dv, levels = c("Disagreement",
                                                "Symbolic-operational\nideological alignment",
                                                "Within-operational\nideological alignment"))

FigureI2 <- ggplot(plotdf_I, aes(x = predictor, y = coef, color = dv)) +
  facet_wrap(~dv, scales = "free_x") +
  geom_pointrange(aes(ymin = coef - 1.96 * se, ymax = coef + 1.96 * se),
                  size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = label), vjust = -0.8, hjust = 0.5,
            size = 3, color = "black") +
  coord_flip() +
  labs(y = "Standardised coefficient", x = "",
       subtitle = "SI Figure I2: Weighted results") +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10), legend.position = "none",
        strip.background = element_rect(fill = "grey90"),
        strip.text = element_text(size = 9))
ggsave("SI_Figure_I2_coefplot_weighted.png",
       plot = FigureI2, height = 5, width = 9)

# Table I1: OLS predicting Disagreement (weighted)
TI1_m1 <- lm(disagreement ~ symbolic_op_align,             data = countrydat_I_scaled)
TI1_m2 <- lm(disagreement ~ symbolic_op_align + HDI + HDI2, data = countrydat_I_scaled)
TI1_m3 <- lm(disagreement ~ within_op_align,               data = countrydat_I_scaled)
TI1_m4 <- lm(disagreement ~ within_op_align + HDI + HDI2,  data = countrydat_I_scaled)

cat("\n=== SI Table I1: OLS Predicting Disagreement (weighted) ===\n")
texreg::screenreg(
  list(TI1_m1, TI1_m2, TI1_m3, TI1_m4),
  digits = 3,
  custom.model.names = c("(1) Symbolic-op",
                         "(2) Symbolic-op + HDI",
                         "(3) Within-op",
                         "(4) Within-op + HDI"),
  custom.coef.names  = c("(Intercept)",
                         "Symbolic-operational ideological alignment",
                         "Within-operational ideological alignment",
                         "HDI", "HDI-squared")
)
cat("Saved SI Section I figures and tables.\n")


################################################################################
# END OF REPLICATION SCRIPT

